Excitonic signatures of ferroelectric order in parallel-stacked MoS2

Interfacial ferroelectricity, prevalent in various parallel-stacked layered materials, allows switching of out-of-plane ferroelectric order by in-plane sliding of adjacent layers. Its resilience against doping potentially enables next-generation storage and logic devices. However, studies have been limited to indirect sensing or visualization of ferroelectricity. For transition metal dichalcogenides, there is little knowledge about the influence of ferroelectric order on their intrinsic valley and excitonic properties. Here, we report direct probing of ferroelectricity in few-layer 3R-MoS2 using reflectance contrast spectroscopy. Contrary to a simple electrostatic perception, layer-hybridized excitons with out-of-plane electric dipole moment remain decoupled from ferroelectric ordering, while intralayer excitons with in-plane dipole orientation are sensitive to it. Ab initio calculations identify stacking-specific interlayer hybridization leading to this asymmetric response. Exploiting this sensitivity, we demonstrate optical readout and control of multi-state polarization with hysteretic switching in a field-effect device. Time-resolved Kerr ellipticity reveals direct correspondence between spin-valley dynamics and stacking order.


S3. Reflectance measurements and spatial filtering
All optical measurements were performed in flow-type cryostats with an accessible temperature range of 4K to room temperature.The samples were kept under vacuum for all the optical measurements.All the experiments were performed in back-scattering geometry.
For white light reflectance measurements, the samples were illuminated with a quartz tungsten halogen lamp.The collimated beam was focused on the sample using an 80X objective.The reflected light from the sample was collected using the same objective and spectrally resolved using a spectrometer and a charge-coupled detector.To obtain the hyperspectral map of the samples, the cryostat with the samples inside, was moved with respect to the fixed optical spot using a computer-controlled XY stage.The data was analyzed with LabView and Origin software.
In the detection path, before the spectrograph, we introduce a home-built Spatial-Filtering module to enhance the spatial resolution.The input side of the spatial filter consists of an aspheric lens mounted on a Z-translator.It focuses the reflected beam onto a pinhole.The pinhole was carefully aligned to the beam path with the help of an XY translator.A fraction of the focused light passes through the pinhole aperture.Another aspheric lens was used to collect the beam from the pinhole and collimate along the detection path.

S4. Time resolved reflectivity and Kerr ellipticity
In the pump-probe setup, two separately tunable pulsed lasers (Toptica: femtoFiberPro) were used, one for the pump and the other for the probe laser beam.Each system emitted with a pulse repetition rate of 80 MHz, a spectral width of 6 meV, a pulse duration of about ∼200 fs.Probe and pump pulses were electronically timesynchronized and amplitude-modulated with different chopping frequencies, adding up to a sum modulation frequency of ∼950 Hz.A mechanical delay line changed the path length of both beams and, thus, the time offset of both pulses.Through an achromatic λ/4 plate, the pump beam was circularly polarized to either left-or right-handed helicity.After reflection the pump beam was filtered out by a long-pass filter while the polarization and signal amplitude of the probe beam were measured.The polarization of the probe beam was analyzed for its ellipticity by a combination of a λ/4 plate, a Wollaston prism, and two photodiodes (ThorLabs PDB210A Si Photodetector).The difference signal of the two diodes was fed into a lock-in amplifier, yielding a TRKE signal.The sum signal of the diodes was obtained by using an external adder and fed into a second lock-in amplifier, yielding the TR∆R signal.Both lock-in amplifiers were fed the same reference frequency given by the modulation sum frequency.With this approach, it was possible to measure the TRKE and TR∆R at the same time, cutting the needed measurement time in half, compared to, as commonly performed, subsequent measurements.The TRKE setup is described in detail elsewhere. 2,3

S5. Density functional theory
To investigate the electronic properties and optical selection rules of the five-layer (5L) MoS 2 with 3R stacking, we performed density functional theory (DFT) calculations using the all-electron full-potential linearized augmented plane-wave (LAPW) method implemented in the Wien2k code. 4 We consider the Perdew-Burke-Ernzerhof 5 exchange-correlation functional with van der Waals interactions included via the D3 correction. 6e wave function expansion into atomic spheres uses orbital quantum numbers up to 10 and the plane-wave cut-off multiplied with the smallest atomic radii is set to 8. Spin-orbit coupling was included fully relativistically for core electrons, while valence electrons were treated within a second-variational procedure 7 with the scalar-relativistic wave functions calculated in an energy window up to 2 Ry.Self-consistency was achieved using a two-dimensional Monkhorst-Pack k-grid with 15×15×1 points and a convergence criteria of 10 −6 e for the charge and 10 −6 Ry for the energy.We used the structural parameters taken from Ref. 8, i. e., the in-plane lattice parameter is 3.191 Å, the thickness of a single MoS 2 layer is 3.116 Å, the interlayer distance between adjacent MoS 2 layers is 3.094 Å, and a vacuum region of 20 Å was used to avoid interaction between the heterostructures' replicas.We explored all the possible combinations of 3R stackings in a 5L system leading to a total of 10 different structures, shown in Fig. S18 (see also Fig. S12).The total energy of these different systems are nearly the same, as shown in Fig. S19, supporting the appearance of different domains in the studied samples.The calculated band structures with spin-orbit coupling along the Γ − K − M line are shown in Fig. S20 with the color-code representing the spin expectation value in the out-of-plane direction.A zoom of the conduction and valence bands in the vicinity of the K point is shown in Fig. S21, revealing that the energy bands that generate A and B excitons show strong spin-valley locking.The wave functions presented in Fig. 2 of the main text were evaluated without spin-orbit coupling.Turning to the selection rules at the K point, we present in Fig. S22 the absolute value of the dipole matrix elements (|p cv | = ℏ m0 | ⟨c, K |⃗ p • α| v, K⟩ | calculated within the LAPW basis set, 9 with α = s + , s − , z encoding the light polarization) as a function of the band labels.
In Fig. S23 we display the oscillator strength, i. e., |p cv | 2 , that enters the absorption spectra.In Figs.S24 and S25 we present |p cv | and |p cv | 2 as a function of the transition energy only for s + polarization (the dominant contribution due to the stron spin-valley locking).Our results reveal that variation of the transition energies is on the order of 10-20 meV for the different stackings.

S6. Exciton binding energies
To evaluate the effect of the asymmetric dielectric surroundings of the 5L MoS 2 systems, we calculated the binding energies for the intralayer A excitons located at each layer.We employed the effective Bethe-Salpeter equation 10,11 formalism, assuming non-interacting parabolic bands for electrons and holes. 12,13 e consider (from Ref. 15), and ε(air) = 1.We consider each MoS 2 layer to have an effective thickness of 6.232 Å, taken as twice the value of the physical thickness of the TMDC. 16The calculated dielectric constants experienced by the intralayer excitons at each MoS 2 layer are show in Fig. S26.The potential at each layer, l, takes the form −1 , with l = 1, . . ., 5. The resulting exciton binding energies are given in Table S1, showing variations of up to 30 meV, a similar energy scale of the optical transitions discussed in Figs.S24 and S25.
For the structures with mixed layer contributions (MX-XM-MX-XM, XM-MX-XM-MX, and XM-XM-MX-MX) we average the energy transitions and leave the intralayer exciton potentials intact.Since these quantities (transition energies and exciton binding energies) operate on the similar energy scales, our conclusions remain valid.Thus, combining the DFT transition energies and dielectric effects within the exciton picture, we present the resulting absorption exciton spectra in Fig. S27.We used the same value of the dipole matrix element for all transitions.DFT underestimates the band gap and therefore our calculations are red-shifted by ∼480 meV in comparison with the experimental values.Nonetheless, the exciton absorption peaks are 10-30 meV apart depending on the particular stacking (ferroelectric domains), in excellent agreement with the experimental findings shown in Fig. 1d of the main text.

S7. Extended data
the average value of effective masses taken from the DFT calculations (m v = 0.515 and m c = 0.415) and the electrostatic potential for each layer is obtained numerically by solving the Poisson equation in k-space, assuming the different regions to have dielectric constants of ε(MoS 2 ) = 15.45 (from Ref. 14), ε(hBN) = 4.5

Figure S1 |
Figure S1 | Multi-polarization states in naturally grown 3R-MoS2.a. Surface potential map copied from Fig.1b of the main text for reference.b.Typical line cut of the lateral surface potential profile across a domain wall separating regions I and II, marked by the red line in panel a.The measured potential variation, ∆ VKP ∼ 127 mVagrees well with earlier reported values in this material.17,18

Figure S2 |
Figure S2 | Inducing ferroelectric domains in naturally grown 3R-MoS2.a. Optical micrograph of a bi-and trilayer 3R-MoS2 flake stamped on Si/SiO2/hBN with deliberate shear perturbation (horizontal force) during the transfer from PDMS. b.Surface potential map within the area enclosed by the black rectangle in a. Densely packed ferroelectric domains can be identified by the difference in contrast.c.A typical line cut of the lateral surface potential profile across domains in the bilayer region, marked by the blue line in panel b.The red circle is the surface potential from the trilayer region.

Figure S3 |
Figure S3 | Reflectance spectra.a. Optical micrograph of Sample 1. b-d.As recorded reflectance spectra at 4K from various spatial locations, as indicated by symbols and lines.

Figure S4 |
Figure S4 | Cavity effect on reflection contrast intensity.Room temperature reflectance contrast spectra of four different samples along with their optical microscopic images.Samples with a turquoise appearance (left panels) shows a pronounced Γ-Q or K-Q and XA peaks.Conversely, sample prepared on thicker hBN, appearing pink (central panels), exhibit no detectable ∆R/R intensity in the lower energy range, and the XA transition is only weakly visible.However, the XB feature is relatively stronger in this sample.In the case of an orange-colored sample (top-encapsulated, right panel), XA and XB transitions are prominently present, while the lower energy feature remains invisible.Here we have used RC = ∆R/R F lake .

Figure S5 |
Figure S5 | Cavity effect on reflection contrast intensity.The same as in Fig. S4 with the exception of RC being equal to ∆R/R Ref .

Figure S6 |
Figure S6 | Comparison of (R Ref − R Flake )/R Flake with (R Ref − R Flake )/R Ref for Sample 1. Integrated intensity map of ∆R/R at the XA spectral region, i.e., from 1906 to 1918 meV along with low-temperature reflectance contrast spectra from various spatial locations marked by symbols of corresponding color in the map (black -5L domain I, red 5L domain II, and blue -7L).We have used in (a-c.)RC = (R Ref − R F lake )/R F lake (as in Fig. 1c and d) and in (d-f.)RC = (R Ref − R F lake )/R Ref .

Figure S7 |
Figure S7 | Reflection contrast intensity of ferroelectric domains in few-layer 3R-MoS2.a-c.Lowtemperature integrated intensity map of ∆R/R of the complete flake at the Γ-Q or K-Q peak, XA, and XB spectral region, respectively.

Figure S8 |
Figure S8 | High resolution reflection contrast intensity of ferroelectric domains in few-layer 3R-MoS2.a-b.Low-temperature integrated intensity map of ∆R/R of the selected region of flake as Fig. 1 in the main text at the low energy tail and Bx spectral region, respectively.This data set differs from Fig. S7 in terms of its spatial and spectral resolution.

Figure S9 |
Figure S9 | First derivative of reflection contrast at room temperature.a. Numerically calculated first derivative of reflectance contrast spectra from Domain I and II.b.Sum intensity map of d dE ∆R R at the XA spectral region i.e. from 1843 to 1863 meV (gray bar in a).The analysis highlights that the ferroelectric domains can be probed at room temperature, opening up an unprecedented opportunity for potential applications.

Figure S10 |
Figure S10 | Reflection contrast imaging of ferroelectric domains in trilayer 3R-MoS2.a. Optical micrograph of a 3R-MoS2 flake on Si/SiO2/hBN.b.Surface potential map of the entire flake.c.Room temperature sum intensity map of d dE ∆R R at the XA spectral region i.e. from 1857 to 1863 meV.d.Low-temperature reflectance contrast map at the XA spectral region i.e. from 1898.2 to 1906.2 meV.

Figure S11 |
Figure S11 | Cavity effect on reflection contrast intensity.a.Another 3R-MoS2 sample on a 95 nm-hBN/285 nm-SiO2/Si substrate.b-d.Low-temperature reflectance contrast spectra from various spatial locations from the sample shown in panel a, plotted in solid line.Due to the chosen thickness of hBN and SiO2 substrate, the reflection contrast at the high energy side increases drastically.The spectra from Fig.1-Sample 1 are included as scatter plots for comparison.

Figure S12 |
Figure S12 | Possible permutations of stacking order in a five-layer flake.We replace the crystallographic depiction with their respective polarization vector for convenience, for instance, an XM stacking configuration with an upwards black colored arrow and MX as a red downwards arrow.Here, we have taken a systematic approach of depicting all the possible combinations by changing the stacking order of one interface at a time, starting from copolarized XM-XM-XM-XM (1) configuration to entirely flipped MX-MX-MX-MX (1 ′ ) ordering.Orange-colored boxes indicate crystallographically non-equivalent and unique stacking orders.The "′" symbol has been used to indicate equivalent yet flipped configurations.The scheme shows all sixteen possible stacking combinations, out of which ten are crystallographically unique, and the rest are flipped versions of a few.For example, XM-MX-MX-XM (8 ′ ) is a 180 o flipped version of MX-XM-XM-MX(8).It is interesting to note that KPFM measurement can distinguish among different rows of this schematic representation but is insensitive to the elements within a given row.For example, XM-XM-XM-XM (1) and XM-XM-XM-MX (2) would have different surface potential; however, all the combinations in the second row, from 2 to 5, would remain indistinguishable in a KPFM measurement.

Figure S13 |
Figure S13 | Schematic of real-space band alignment for different stacking, a purely electrostatic perspective.Schematics of layer-projected KVB and KCB band-edges for crystallographically unique stacking configurations, discussed in Fig. S12.The dashed lines act as a guide to the eye to highlight valence-band degeneracies, which can occur in partially polarized stacking.To first order, the band-edge variation arises due to the rigid shift in layer-projected K valleys introduced by the ferroelectricity-induced electrostatic considerations.

Figure S14 |
Figure S14 | Characterization of Sample 4. a. Optical micrograph of a top-encapsulated 3L-3R-MoS2 flake on Si/SiO2/hBN.Topographical steps and edges of the MoS2 flake have been marked by white lines.We have used blue lines to mark the graphite flake used as a source/drain contact.On the right panel we draw the device schematic from a lateral perspective.b.Surface potential map at different gate voltages within the area enclosed by the black rectangle in a.Two domains can be identified by the difference in contrast, which is particularly strong for Vg= -8 V.

Figure S15 |
Figure S15 | Comparison of (R Ref − R Flake )/R Flake with (R Ref − R Flake )/R Ref for Sample 4. a -e Replica of Fig. 3 i.e.RC maps, RC intensity as a function of gate voltage, RC spectra, and its hysteretic response where we have used RC=(R Ref − R F lake )/R F lake .(f -j) the same with RC defined as (R Ref − R F lake )/R Ref ,

Figure S17 |
Figure S17 | Integrated RC intensity (≡ line cuts from panel b and c of Fig. 3) and standard deviation profile at the A and B excitonic region.Integrated intensity = xmax x=x min ∆R/R(x); Standard deviation = xmax x=x min ∆R/R(x)(x−x 0 ) 2 xmax x=x min ∆R/R(x)

Figure S20 |Figure S21 |
Figure S20 | Band structure along the Γ − K − M line, color-coded according to the expectation value of spin along Z, < sz >.

Figure S22 |
Figure S22 | Absolute value of the dipole matrix elements between the valence bands with spin-up (different columns) to all conduction bands (labeled by index).

Figure
Figure S28 | a-b.False color map of transient ellipticity as a function of probe energy from domain I and II, respectively.c-d.Transient reflection at selected energies, marked by colored arrows in the top panels (copied from Fig. 4 of the main text for completeness).